#-*-coding:utf-8-*-


from csv import Dialect, DictReader
from distutils.command.build_scripts import first_line_re
import encodings
from operator import truediv
from pdb import pm
from tkinter.tix import DirSelectBox
from turtle import backward
import numpy as np
#from pictureplot import pictureplotting
import random
import matplotlib.pyplot as plt
import platform
import os
from plot import plottree

import global_var


def all(filename):
    """读取filename所指.dat文件进行潮流计算和承载力评估
    返回包含各节点计算结果的字典datadict
    """
    def inputdata(path):
        """得到文件所有数据行的内容
        
        参数
        -----------------
        path : string
            所选.dat文件在当前用户端的路径

        返回值
        -----------------
        data : list
            data的元素是.dat文件每一数据行的string
        """
        data = []
        with open(path, 'r', encoding='gbk') as f:
            for line in f.readlines():
                if line.startswith('B') or line.startswith('L') or line.startswith('T'):
                    data.append(line)
        return data

    def get_bcard(data):
        """获得节点数据卡的各项数据

        参数
        -----------------
        data : list
            包含所有数据行

        返回值
        -----------------
        bdata : 2D list
            元素为一个B节点数据形成的list
        """
        bdata = []
        for line in data:
            mline = []
            if line.startswith('B'):
                line = line.encode("gbk")

                mline.append(line[1:2].strip().decode(encoding='gbk',errors='ignore'))  # 卡片子型
                mline.append(line[2:3].strip().decode(encoding='gbk',errors='ignore'))  # 修改码
                mline.append(line[3:6].strip().decode(encoding='gbk',errors='ignore'))  #所有者代码
                mline.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))  # 节点名称
                mline.append(line[14:18].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压(kV)
                mline.append(line[18:20].strip().decode(encoding='gbk',errors='ignore'))  # 节点所在分区名称
                mline.append(line[20:25].strip().decode(encoding='gbk',errors='ignore'))  # 有功恒定负荷(MW) 
                mline.append(line[25:30].strip().decode(encoding='gbk',errors='ignore'))  # 无功恒定负荷(MVar) 
                mline.append(line[30:34].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压下的有功并联导纳负荷(MW) 
                mline.append(line[34:38].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压下的无功并联导纳负荷(MVar)
                mline.append(line[38:42].strip().decode(encoding='gbk',errors='ignore'))  # 最大有功出力(MW)
                mline.append(line[42:47].strip().decode(encoding='gbk',errors='ignore'))  # 实际有功出力(MW) 
                mline.append(line[47:52].strip().decode(encoding='gbk',errors='ignore'))  # （B）所安排的无功出力值(MVar)/（BQ、BS）无功出力最大值(MVar)
                mline.append(line[52:57].strip().decode(encoding='gbk',errors='ignore'))  # 无功出力最小值(MVar)
                mline.append(line[57:61].strip().decode(encoding='gbk',errors='ignore'))  # 所安排的电压值/Vmax(标幺值) 
                mline.append(line[61:65].strip().decode(encoding='gbk',errors='ignore'))  # （B、BQ）所安排的Vmin(标幺值) /（BS）角度值

                # print(mline)
                bdata.append(mline)
        # bdata_array = np.array(bdata)
        # print(bdata_array)
        return bdata

    def get_lcard(data):
        """获得对称线路数据卡的各项数据

        参数
        -----------------
        data : list
            包含所有数据行

        返回值
        -----------------
        ldata_array : 2D numpy.array
            一行包含一个L线路的数据
        """
        ldata = []
        for line in data:
            mline = []
            if line.startswith('L'):
                line = line.encode("gbk")

                mline.append('')
                mline.append(line[2:3].strip().decode(encoding='gbk',errors='ignore'))  # 修改码
                mline.append(line[3:6].strip().decode(encoding='gbk',errors='ignore'))  #所有者代码
                mline.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))  # 节点名1
                mline.append(line[14:18].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压(kV)
                mline.append(line[18:19].strip().decode(encoding='gbk',errors='ignore'))  # 区域联络线测点标记
                mline.append(line[19:27].strip().decode(encoding='gbk',errors='ignore'))  # 节点名2
                mline.append(line[27:31].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压(kV) 
                mline.append(line[31:32].strip().decode(encoding='gbk',errors='ignore'))  # 回路号
                mline.append(line[33:37].strip().decode(encoding='gbk',errors='ignore'))  # 线路额定电流值(A) 
                mline.append(line[37:38].strip().decode(encoding='gbk',errors='ignore'))  # 并联线路数目
                mline.append(line[38:44].strip().decode(encoding='gbk',errors='ignore'))  # 系统基准电压和基准容量下的电阻标幺值
                mline.append(line[44:50].strip().decode(encoding='gbk',errors='ignore'))  # 系统基准电压和基准容量下的电抗标幺值
                mline.append(line[50:56].strip().decode(encoding='gbk',errors='ignore'))  # 线路对地导纳标幺值G/2
                mline.append(line[56:62].strip().decode(encoding='gbk',errors='ignore'))  # 线路对地导纳标幺值B/2
                mline.append(line[62:66].strip().decode(encoding='gbk',errors='ignore'))  # 线路长度
                mline.append(line[66:74].strip().decode(encoding='gbk',errors='ignore'))  # 线路说明数据
                mline.append(line[74:75].strip().decode(encoding='gbk',errors='ignore'))  # 投运日期：月
                mline.append(line[75:77].strip().decode(encoding='gbk',errors='ignore'))  # 投运日期：年
                mline.append(line[77:80].strip().decode(encoding='gbk',errors='ignore'))  # 停运日期

                # print(mline)
                ldata.append(mline)
        ldata_array = np.array(ldata)
        # print(ldata_array)
        return ldata_array

    def get_tcard(data):
        """获得变压器/移相器数据卡的各项数据

        参数
        -----------------
        data : list
            包含所有数据行

        返回值
        -----------------
        tdata_array : 2D numpy.array
            一行包含一个T变压器/TP移相器的数据
        """
        tdata = []
        for line in data:
            mline = []
            if line.startswith('T'):
                line = line.encode("gbk")

                mline.append('T') # 卡片类型
                mline.append(line[1:2].strip().decode(encoding='gbk',errors='ignore'))  # 卡片子型
                mline.append(line[2:3].strip().decode(encoding='gbk',errors='ignore'))  # 修改码
                mline.append(line[3:6].strip().decode(encoding='gbk',errors='ignore'))  #所有者代码
                mline.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))  # 节点名1
                mline.append(line[14:18].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压(kV)
                mline.append(line[18:19].strip().decode(encoding='gbk',errors='ignore'))  # 区域联络线测点标记
                mline.append(line[19:27].strip().decode(encoding='gbk',errors='ignore'))  # 节点名2
                mline.append(line[27:31].strip().decode(encoding='gbk',errors='ignore'))  # 基准电压(kV) 
                mline.append(line[31:32].strip().decode(encoding='gbk',errors='ignore'))  # 回路号
                mline.append(line[33:37].strip().decode(encoding='gbk',errors='ignore'))  # 变压器额定容量 
                mline.append(line[37:38].strip().decode(encoding='gbk',errors='ignore'))  # 并联线路数目
                mline.append(line[38:44].strip().decode(encoding='gbk',errors='ignore'))  # 由铜损引起的等效电阻标幺值
                mline.append(line[44:50].strip().decode(encoding='gbk',errors='ignore'))  # 漏抗标幺值
                mline.append(line[50:56].strip().decode(encoding='gbk',errors='ignore'))  # 由铁损引起的等效电导标幺值
                mline.append(line[56:62].strip().decode(encoding='gbk',errors='ignore'))  # 激磁电纳标幺值
                mline.append(line[62:67].strip().decode(encoding='gbk',errors='ignore'))  # （T）节点1的分接头位置(kV)/（TP）移相器的固定移相角(度)
                mline.append(line[67:72].strip().decode(encoding='gbk',errors='ignore'))  # （T）节点2的分接头位置(kV)
                mline.append(line[74:75].strip().decode(encoding='gbk',errors='ignore'))  # 投运日期：月
                mline.append(line[75:77].strip().decode(encoding='gbk',errors='ignore'))  # 投运日期：年
                mline.append(line[77:80].strip().decode(encoding='gbk',errors='ignore'))  # 停运日期

                # print(mline)
                tdata.append(mline)
        tdata_array = np.array(tdata)
        # print(tdata_array)
        return tdata_array

    def get_mx():
        '''
        获取发电机次暂态电抗

        参数
        -----------------
        nodename : list
            包含所有节点名称的全局变量

        返回值
        -----------------
        mxdata : 2D list
            节点的数据行编号和暂态电抗
        '''
        mxdata = []
        for i, name in enumerate(nodename):
            if name == '4':
                xd = [i + 1, 0.102]
                mxdata.append(xd)
            elif name == '50':
                xd = [i + 1, 0.072]
                mxdata.append(xd)
            elif name == '100':
                xd = [i + 1, 0.01]
                mxdata.append(xd)
        # print(mxdata)
        return mxdata

    def get_nodedata(data):
        """获得各节点数据

        参数
        -----------------
        data : list
            包含所有数据行

        返回值
        -----------------
        nodedata_array : 2D numpy.array
            一行包含一个节点的数据: bus_i type Pd Qd Pg Qg Gs Bs Vm Va
        nodedata : 2D list
            用于传入get_linedata()内部继续使用
        """
        nodedata = []
        for line in data:
            mline = []
            if line.startswith('B'):
                line = line.encode("gbk")

                mline.append(line[0:2].strip().decode(encoding='gbk',errors='ignore'))  # card id, syb type 
                mline.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))  # name 
                mline.append(line[14:18].strip().decode(encoding='gbk',errors='ignore'))  # kV 
                mline.append(line[20:25].strip().decode(encoding='gbk',errors='ignore'))  # PL(MW) 
                mline.append(line[25:30].strip().decode(encoding='gbk',errors='ignore'))  # QL(MVar) 
                mline.append(line[30:34].strip().decode(encoding='gbk',errors='ignore'))  # Pshunt(MW) 
                mline.append(line[34:38].strip().decode(encoding='gbk',errors='ignore'))  # Qshunt(MVar) 
                mline.append(line[42:47].strip().decode(encoding='gbk',errors='ignore'))  # PGen(MW) 
                mline.append(line[57:61].strip().decode(encoding='gbk',errors='ignore'))  # Vmax(V) 

                nodename.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))
                nodeulv.append(line[14:18].strip().decode(encoding='gbk',errors='ignore'))

                # print(mline)
                nodedata.append(mline)
        # print(nodename)
        # print(nodeulv)
        nodedata_array = np.zeros([len(nodedata), 10])
        for i in range(len(nodedata)):
            nodedata_array[i, 0] = int(i + 1)  # bus_i: 从1开始按顺序编号节点

            # B : PQ节点 -> 1
            # BQ : PV节点 -> 2
            # BS : 平衡节点 -> 3
            if nodedata[i][0] == 'BQ':
                nodedata_array[i, 1] = 2
            elif nodedata[i][0] == 'BS':
                nodedata_array[i, 1] = 3
            else:
                nodedata_array[i, 1] = 1

            # Pd: 节点负荷有功功率
            if nodedata[i][3] == '':
                nodedata_array[i, 2] = 0
            else:
                nodedata_array[i, 2] = float(nodedata[i][3]) / 100
            # Qd: 节点负荷无功功率
            if nodedata[i][4] == '':
                nodedata_array[i, 3] = 0
            else:
                nodedata_array[i, 3] = float(nodedata[i][4]) / 100
            # Pg: 节点发电机输出有功功率
            if nodedata[i][5] == '':
                nodedata_array[i, 4] = 0
            else:
                nodedata_array[i, 4] = float(nodedata[i][5]) / 100
            # PV_Pg: PV节点的节点发电机输出有功功率
            if nodedata[i][0] == 'BQ':
                nodedata_array[i, 4] = float(nodedata[i][-2]) / 100
            # Qg: 节点发电机输出无功功率
            if nodedata[i][6] == '':
                nodedata_array[i, 5] = 0
            else:
                nodedata_array[i, 5] = float(nodedata[i][6]) / 100
            # Vm: 电压幅值
            if nodedata[i][-1] == '':
                nodedata_array[i, 8] = 1.
            else:
                nodedata_array[i, 8] = float(nodedata[i][-1])
            nodedata_array[i, 6] = 0  # Gs: 节点对地电导
            nodedata_array[i, 7] = 0  # Bs: 节点对地导纳
            nodedata_array[i, 9] = 0  # Va: 电压相角
        # nodedata_array = np.array(nodedata_array)
        # print(nodedata_array)
        return nodedata_array,nodedata

    def get_linedata(data, nodedata_origin):
        """获得各线路及变压器的数据

        参数
        -----------------
        data : list
            包含所有数据行
        nodedata_origin : 2D list
            来自get_nodedata()对数据行的初步处理

        返回值
        -----------------
        linedata_array : 2D numpy.array
            一行包含一条线路或一台变压器的数据: fbus  tbus   r    x    b    ratio
        """
        k = 1  
        nodedata = []
        linedata = []
        for mline in nodedata_origin:
            mline.append(k)  # k: 从1开始按顺序编号线路
            k += 1
            nodedata.append(mline)
        for line in data:
            nline = []  # 线路数据
            oline = []  # 变压器数据
            line = line.encode("gbk")

            #处理线路数据
            if line.decode(encoding='gbk',errors='ignore').startswith('L'):
                # fbus: 起始节点
                for mdata in nodedata:
                    # print(mdata)
                    if mdata[1] == line[6:14].strip().decode(encoding='gbk',errors='ignore'):
                        nline.append(int(mdata[-1]))
                # tbus: 终止节点
                for mdata in nodedata:
                    if mdata[1] == line[19:27].strip().decode(encoding='gbk',errors='ignore'):
                        nline.append(mdata[-1])
                        # r: 支路电阻
                        if line[38:44].strip().decode(encoding='gbk',errors='ignore') != '':
                            nline.append(float(line[38:44].strip().decode(encoding='gbk',errors='ignore')))
                        else:
                            nline.append(0.0)
                        # x: 支路电抗
                        if line[44:50].strip().decode(encoding='gbk',errors='ignore') != '':
                            nline.append(float(line[44:50].strip().decode(encoding='gbk',errors='ignore')))
                        else:
                            nline.append(0.0)
                        # b: 节点对地导纳
                        if line[56:62].strip().decode(encoding='gbk',errors='ignore') != '':
                            nline.append(float(line[56:62].strip().decode(encoding='gbk',errors='ignore')))
                        else:
                            nline.append(0.0)
                        nline.append(0)  # 线路没有变比属性
                # print(nline)
                linedata.append(nline)
            
            #处理变压器数据
            if line.decode(encoding='gbk',errors='ignore').startswith('T'):
                # fbus: 起始节点
                for mdata in nodedata:
                    if mdata[1] == line[6:14].strip().decode(encoding='gbk',errors='ignore'):
                        oline.append(int(mdata[-1]))
                # tbus: 终止节点
                for mdata in nodedata:
                    if mdata[1] == line[19:27].strip().decode(encoding='gbk',errors='ignore'):
                        oline.append(mdata[-1])
                        oline.append(float(line[38:44].strip().decode(encoding='gbk',errors='ignore')))  # r: 支路电阻
                        oline.append(float(line[44:50].strip().decode(encoding='gbk',errors='ignore')))  # x: 支路电抗
                        oline.append(0.0)  # 变压器没有对地导纳属性
                        oline.append(1.05)  # ratio: 变比
                # print(oline)
                linedata.append(oline)

        # print(liendata)
        linedata_array = np.array(linedata)
        return linedata_array
    
    def get_netdata(nodedata):
        """以bus_i为标签分类节点并计算输入功率

        参数
        -----------------
        nodedata : 2D numpy.array
            一行包含一个节点的数据: bus_i type Pd Qd Pg Qg Gs Bs Vm Va

        返回值
        -----------------
        PQnode : 2D list
            PQ节点的bus_i
        PVnode : 2D list
            PV节点的bus_i
        slacknode : 2D list
            平衡节点的bus_i
        P_real : list
            按bus_i编号顺序列出输入有功功率
        Q_real : list
            按bus_i编号顺序列出输入无功功率
        """
        PQnode = nodedata[np.where(nodedata[:, 1] == 1), 0]  # PQ节点
        PVnode = nodedata[np.where(nodedata[:, 1] == 2), 0]  # PV节点
        slacknode = nodedata[np.where(nodedata[:, 1] == 3), 0]  # 平衡节点
        P_real = -nodedata[:, 2] + nodedata[:, 4]  # 节点输入有功功率
        Q_real = -nodedata[:, 3] + nodedata[:, 5]  # 节点输入无功功率
        # print(type(PQnode))
        # print(PVnode)
        # print(slacknode)
        # print(P_real)
        # print(Q_real)
        return (PQnode, PVnode, slacknode, P_real, Q_real)

    def get_y(nodedata, linedata, mdata):
        """获得导纳矩阵元素

        参数
        -----------------
        nodedata : 2D numpy.array
            一行包含一个节点的数据: bus_i type Pd Qd Pg Qg Gs Bs Vm Va
        linedata : 2D numpy.array
            一行包含一条线路或一台变压器的数据: fbus  tbus   r    x    b    ratio
        mdata : 2D list
            节点的数据行编号和暂态电抗

        返回值
        -----------------
        Y : 2D numpy.array
            导纳矩阵元素
        """
        numnode = nodedata.shape[0]
        numline = linedata.shape[0]
        Y = np.zeros([numnode, numnode]) + np.zeros([numnode, numnode]) * 1j
        G0 = nodedata[:, 6]  # 节点对地电导
        B0 = nodedata[:, 7]  # 节点对地电纳
        for i in range(numline):
            node1 = int(linedata[i][0] - 1)
            node2 = int(linedata[i][1] - 1)
            # print(node1,node2)
            R = linedata[i][2]
            X = linedata[i][3]

            if linedata[i][5] == 0:  # 普通线路，无变压器
                B_2 = linedata[i][4]
                Y[node1, node1] = Y[node1, node1] + B_2 * 1j + 1 / (R + 1j * X)
                Y[node2, node2] = Y[node2, node2] + B_2 * 1j + 1 / (R + 1j * X)
                Y[node1, node2] = Y[node1, node2] - 1 / (R + 1j * X)
                Y[node2, node1] = Y[node2, node1] - 1 / (R + 1j * X)
            else:  # 有变压器支路
                K = linedata[i][5]
                YT = 1 / (R + 1j * X)
                Y[node1, node1] = Y[node1, node1] + (K - 1) / K * YT + YT / K
                Y[node2, node2] = Y[node2, node2] + (1 - K) / K ** 2 * YT + YT / K
                Y[node1, node2] = Y[node1, node2] - 1 / K * YT
                Y[node2, node1] = Y[node2, node1] - 1 / K * YT

        for i in range(numnode):
            node = int(nodedata[i, 0] - 1)  # 第一列为节点编号
            Y[node, node] = Y[node, node] + G0[i] + 1j * B0[i]
        for data in mdata:
            Y[data[0], data[0]] = Y[data[0], data[0]] + 1 / data[1]
        # print(Y)
        return (Y)

    def polar_u(U_amplitude, U_angle):
        """获得初始电压和相角

        参数
        -----------------
        U_amplitude : 2D numpy.array
            电压幅值
        U_angle : 2D numpy.array
            电压相角

        返回值
        -----------------
        U : numpy.array
            初始电压
        angle : numpy.array
            初始相角
        """
        U_complex = U_amplitude * np.cos(U_angle) + U_amplitude * np.sin(U_angle) * 1j
        U = abs(U_complex)
        angle = np.angle(U_complex)
        # print(type(U))
        return (U, angle)

    # ---------------直角坐标形式--------------
    def RectanU(U_Amplitude, U_Angle):
        U_e = U_Amplitude * np.cos(U_Angle)
        U_f = U_Amplitude * np.sin(U_Angle)
        return (U_e, U_f)

    def polar_nr(U, angle, Y, PQnode, PVnode, slacknode, P_real, Q_real, tol):
        """进行潮流计算分析

        参数
        -----------------
        U : numpy.array
            初始电压
        angle : numpy.array
            初始相角
        Y : 2D numpy.array
            导纳矩阵元素
        PQnode : 2D list
            PQ节点的bus_i
        PVnode : 2D list
            PV节点的bus_i
        slacknode : 2D list
            平衡节点的bus_i
        P_real : list
            按bus_i编号顺序列出输入有功功率
        Q_real : list
            按bus_i编号顺序列出输入无功功率
        tol : float
            误差限度

        返回值
        -----------------
        U : numpy.array
            该时刻所有节点的电压
        angle : numpy.array
            该时刻所有节点的相角
        maxerror : float
            最大误差
        P : 2D numpy.array
            该时刻所有节点的有功功率
        Q : 2D numpy.array
            该时刻所有节点的无功功率
        """
        P_iter = 0  # 为形成雅可比矩阵
        Q_iter = 0  # 为形成雅可比矩阵
        PQnode = PQnode - 1
        PVnode = PVnode - 1
        slacknode = slacknode - 1
        numnode = Y.shape[0]  # 节点数目
        numPQ = max(PQnode.shape)  # PQ节点数目
        G = Y.real
        B = Y.imag
        P = np.zeros([numnode, 1])
        Q = np.zeros([numnode, 1])
        deltaP = np.zeros([numnode - 1, 1])
        deltaQ = np.zeros([numPQ, 1])
        for i in range(numnode):  # 求解功率不平衡量
            P[i] = U[i] * np.sum(U * (G[i, :] * np.cos(angle[i] 
                                 - angle) + B[i, :] * np.sin(angle[i] - angle)))
            Q[i] = U[i] * np.sum(U * (G[i, :] * np.sin(angle[i] 
                                 - angle) - B[i, :] * np.cos(angle[i] - angle)))
            if i != slacknode:  # 不是平衡节点
                deltaP[P_iter] = P_real[i] - P[i]  # numPQ+numPV
                if i in PQnode:  # PQ节点
                    deltaQ[Q_iter] = Q_real[i] - Q[i]  # NumPQ
                    Q_iter = Q_iter + 1
                P_iter = P_iter + 1
        # print(P)
        # print(Q)
        deltaPQ = np.vstack([deltaP, deltaQ])  # 功率不平衡量

        maxerror = np.max(np.abs(deltaPQ))
        # print(maxerror)
        if maxerror < tol:
            return (U, angle, maxerror, Pi, Qi)

        ## H and N
        HN_iter = -1  # 初始化雅可比矩阵
        H = np.zeros([numnode - 1, numnode - 1])
        N = np.zeros([numnode - 1, numPQ])
        # H and N
        for i in range(numnode):
            if i != slacknode:  # PQ或PV节点
                H_iter_y = -1
                N_iter_y = -1
                HN_iter = HN_iter + 1  # 记录H和N的行数
                for j in range(numnode):
                    if j != slacknode:
                        H_iter_y = H_iter_y + 1  # 记录H列数
                        if i != j:  # 非平衡节点计算H矩阵
                            angleij = angle[i] - angle[j]
                            H[HN_iter, H_iter_y] = -U[i] * U[j] * (
                                G[i, j] * np.sin(angleij) - B[i, j] * np.cos(angleij))
                        else:
                            H[HN_iter, H_iter_y] = Q[i] + U[i] ** 2 * B[i, i]
                        if j in PQnode:
                            N_iter_y = N_iter_y + 1  # 记录N的列数
                            if i != j:
                                angleij = angle[i] - angle[j]
                                N[HN_iter, N_iter_y] = -U[i] * U[j] * (
                                    G[i, j] * np.cos(angleij) + B[i, j] * np.sin(angleij))
                            else:
                                N[HN_iter, N_iter_y] = -P[i] - G[i, i] * U[i] ** 2

        ## J and L
        JL_iter = -1  # 初始化雅可比矩阵
        J = np.zeros([numPQ, numnode - 1])
        L = np.zeros([numPQ, numPQ])
        for i in range(numnode):
            if i in PQnode:  # PQ节点
                JL_iter = JL_iter + 1  # J和L的行数
                J_iter_y = -1
                L_iter_y = -1
                for j in range(numnode):
                    if j != slacknode:  # 非平衡节点
                        J_iter_y = J_iter_y + 1
                        if i != j:
                            angleij = angle[i] - angle[j]
                            J[JL_iter, J_iter_y] = U[i] * U[j] * (
                                G[i, j] * np.cos(angleij) + B[i, j] * np.sin(angleij))
                        else:
                            J[JL_iter, J_iter_y] = -P[i] + G[i, i] * U[i] ** 2
                        if j in PQnode:  # PQ节点
                            L_iter_y = L_iter_y + 1
                            if i != j:
                                angleij = angle[i] - angle[j]
                                L[JL_iter, L_iter_y] = -U[i] * U[j] * (
                                    G[i, j] * np.sin(angleij) - B[i, j] * np.cos(angleij))
                            else:
                                L[JL_iter, L_iter_y] = -Q[i] + B[i, i] * U[i] ** 2

        # 修正
        jaccobi = np.vstack([np.hstack([H, N]), np.hstack([J, L])])
        delta = np.linalg.solve(jaccobi, deltaPQ)
        deltaangle = delta[0:numnode - 1]  # 注意下标
        deltaU_U = delta[numnode - 1:]
        DA_iter = -1
        U_U_iter = -1
        for i in range(numnode):
            if i != slacknode:
                DA_iter = DA_iter + 1
                angle[i] = angle[i] - deltaangle[DA_iter]
                if i in PQnode:
                    U_U_iter = U_U_iter + 1
                    U[i] = U[i] - U[i] * deltaU_U[U_U_iter]

        return (U, angle, maxerror, P, Q)

    cnode = []  # 评估节点

    def getdict(data):
        """确定各节点通过线路或变压器的连接关系

        参数
        -----------------
        data : list
            包含所有数据行

        返回值
        -----------------
        root : string
            当前系统的最高级节点
        mdict : dictionary
            各节点(key)的下级节点(value:list)
        """
        for i in range(len(nodename)):
            if float(nodeulv[i]) > 10:
                cnode.append(nodename[i])
        mdict = {}
        dictdata = []
         # slcaknode = '100'
        for line in data:
            l = []
            if line.startswith('T'):
                line = line.encode("gbk")
                l.append(line[19:27].strip().decode(encoding='gbk',errors='ignore'))  #变压器一次侧
                l.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))  #变压器二次侧
                dictdata.append(l)
        for line in data:
            l = []
            if line.startswith('L'):
                line = line.encode("gbk")
                l.append(line[6:14].strip().decode(encoding='gbk',errors='ignore'))  #线路输电侧
                l.append(line[19:27].strip().decode(encoding='gbk',errors='ignore'))  #线路用电侧
                dictdata.append(l)
        root = ''
        for n in nodename:
            isroot = []
            for d in dictdata:
                if n == d[1]:
                    isroot.append(False)
                else:
                    isroot.append(True)
            if False not in isroot:
                root = n
                break
        # print(root)
        for c in cnode:
            l = []
            for d in dictdata:
                if d[0] == c and d[1] in cnode:
                    l.append(d[1])
            mdict[c] = l
        #print(mdict)
        return (root, mdict)

    #输出方式名和项目名
    def output_name(path):
        """输出方式名和项目名

        参数
        -----------------
        path : string
            所选.dat文件在当前用户端的路径

        返回值
        -----------------
        caseid : string
            方式名
        project : string
            项目名
        """
        with open(path,'r',encoding='gbk',errors='ignore') as file:
             for line in file.readlines():
                 if line.startswith('('):
                    line=line.strip('(')
                    line=line.strip()
                    line=line.strip(')')
                    list=line.split(',')
                    break
        # print(line)
        # print(list)
        caseid=list[1][7:]
        project=list[2][8:]
        # print(type(caseid))
        # print(project)
        return (caseid,project)

    iter = 0
    maxiter = 100
    tol = 10e-10
    nodename = []
    nodeulv = []
    nodedata, nodedata_list = get_nodedata(inputdata(filename))
    rootnode, nodedict = getdict(inputdata(filename))
    # print(rootnode)
    # print(nodedict)
    bcarddata_list = get_bcard(inputdata(filename)) #B数据卡不带评估数据
    lcarddata = get_lcard(inputdata(filename))  #L数据卡
    tcarddata = get_tcard(inputdata(filename))  #L数据卡
    # fdict = {}
    #pictureplotting(rootnode, nodedict)  #输出线状图
    U, angle = polar_u(nodedata[:, 8], nodedata[:, 9])
    # print(nodename)
    PQnode, PVnode, slacknode, P_real, Q_real = get_netdata(nodedata)
    Y = get_y(nodedata, get_linedata(inputdata(filename), nodedata_list), get_mx())
    # print(Y)
    
    while True:
        iter = iter + 1
        # print(U)
        # print(angle)
        U, angle, maxerror, Pi, Qi = polar_nr(U, angle, Y, PQnode, PVnode, slacknode, P_real, Q_real, tol)
        if iter > maxiter or maxerror < tol:
            break  # 结束循环，输出结果
    # print(U)
    # print(Angle*180/3.14159)

    # 将功率转换为标幺值
    def calhot(pd, pl, se=1):
        """反向负载率
        
        参数
        -----------------
        pd : float
            电源出力
        pl : float
            同时刻等效用电负荷
        se : float
            变压器或线路实际容量

        返回值
        -----------------
        lamuda : float
            反向负载率
        """
        lamuda = (pd - pl) / se
        return lamuda

    def calmaxcapicity(lamuda, se=1, kr=0.8):
        """节点的最大容量
        
        参数
        -----------------
        lamuda : float
            评估周期内反向负载率的最大值
        se : float
            变压器或线路实际容量
        kr : float
            设备运行裕度系数

        返回值
        -----------------
        pm : float
            最大容量
        """
        pm = (1 - lamuda) * se * kr
        return pm

    # 热稳定评估主要针对变压器
    linenode = get_linedata(inputdata(filename), nodedata_list)
    #linenode = GetLineData(inputdata(r'C:\Users\19926\Desktop\DG\IEEEdist.DAT'))
 
    # 开始热稳定校核
    res1 = []

    cap = []  # 承载力
    capnode = []  # 分布式电源添加点
    Pi[2] = 0.038
    Pi[-2] = 0
    # print(Pi)
    for node in nodedata:
        if float(nodeulv[int(node[0] - 1)]) > 10:
            hn = node[0]
            capnode.append(hn)

            # print(node)
            # print(hn)
            # print(ln)
            hp = Pi[int(hn) - 1, 0]
            # print(hp)

            lamuda = calhot(0, hp)
            res1.append(lamuda)
            # print(lamuda)
            cap.append(calmaxcapicity(lamuda))
            # print(calmaxcapicity(lamuda))

    # 短路电流校核
    # 10kv限值20ka
    def decideI(ixz, im):
        """新能源接入前系统各级母线节点短路电流不超过相应断路器开断电流限值
        
        参数
        -----------------
        ixz : float
            母线节点短路电流
        im : float
            断路器开断电流限值
        
        返回值
        -----------------
        b : bool
            新能源接入前系统各级母线节点短路电流是否不超过相应断路器开断电流限值
        """
        b = ixz < im
        return b

    def decideaddI(ixz, im, pm, un):
        """新能源接入后系统各级母线节点短路电流不超过相应断路器开断电流限值
        
        参数
        -----------------
        ixz : float
            系统最大运行方式下母线短路电流
        im : float
            断路器开断电流限值
        pm : float
            待校核新能源容量
        un : float
            各级母线额定电压

        返回值
        -----------------
        a : bool
            新能源接入后系统各级母线节点短路电流是否不超过相应断路器开断电流限值
        """
        a = abs(ixz+1.5*pm/np.sqrt(3)/un)
        return a

    # 节点短路电流计算
    """
    短路电流
    220  40    100/230    100/13.8
    110  31.5
    66   31.5
    35  25
    10 20
    """
    Z = np.linalg.inv(Y)

    # 100/230
    def getIxz(vf, zff):
        """按母线最大方式短路容量确定系统最大运行方式下母线短路电流
        
        参数
        -----------------
        un : float
            各级母线额定电压
        zff : float
            某节点在阻抗矩阵的自阻抗

        返回值
        -----------------
        ixz : float
            系统最大运行方式下母线短路电流
        """
        ixz = vf / zff
        return ixz

    # 电流限制转化为标幺值
    mim = [88, 34.65, 20.79, 8.75, 2]

    # 电压偏差校核
    def checkU(un, r, x, p, q, q_, du):
        """待校核容量的分布式光伏接入后产生的电压偏差
        
        参数
        -----------------
        un : float
            该区域内额定电压
        r : float
            节点的自阻抗的电阻分量
        x : float
            节点的自阻抗的电抗分量
        p : float
            分布式光伏的有功功率增量
        q : float
            分布式光伏的无功功率增量正最大值
        q_ : float
            分布式光伏的无功功率增量负最大值
        du : float
            根据各级母线电压和国家标准《电能质量供电电压偏差》估算出的新能源接入导致的电压偏差允许值

        返回值
        -----------------
        checku : bool
            电压偏差校核结果
        """
        Du = abs((r * p + x * q) / un / un)
        Du_ = abs((r * p + x * q_) / un / un)
        # print(Du)
        checku = abs(Du) < abs(du) and abs(Du_) < abs(du)
        return checku

    # 求解节点du
    # print(Z)

    # print(cap)
    # print(capnode)
    # 未带入可新增容量前的短路校核和电压偏差校核
    preires = []
    preures = []
    for i in range(len(res1)):
        preires.append(0)
        preures.append(0)
    for node in nodedata:
        nodenum = int(node[0] - 1)
        for i, cnode in enumerate(capnode):
            if cnode == node[0]:
                if nodeulv[nodenum] == '69.0':
                    ixz = getIxz(node[-2], Z[nodenum, nodenum])
                    preires[i] = decideI(ixz, mim[2])
                if nodeulv[nodenum] == '13.8':
                    ixz = getIxz(node[-2], Z[nodenum, nodenum])
                    preires[i] = decideI(ixz, mim[-1])

    for node in nodedata:
        nodenum = int(node[0] - 1)
        # print(node)
        for i, cnode in enumerate(capnode):
            if cnode == node[0]:
                r = abs(np.real(Z[nodenum, nodenum]))
                x = abs(np.imag(Z[nodenum, nodenum]))
                un = 1
                du = 0.1
                preures[i] = checkU(un, r, x, 0, 0, 0, du)

    # print(cap)
    # print(resu)
    # print(capnode)
    # 初始化参数
    isok = True
    res2 = []
    res3 = []
    resu = []
    for i in range(len(res1)):
        res2.append(0)
        res3.append(0)
        resu.append(0)

    # 主要流程循环
    # 可以尝试改进可新增容量递减方式
    while isok:
        isok = False
        for node in nodedata:
            nodenum = int(node[0] - 1)
            for i, cnode in enumerate(capnode):
                if cnode == node[0]:
                    if nodeulv[nodenum] == '69.0':
                        ixz = getIxz(node[-2], Z[nodenum, nodenum])
                        res2[i] = decideaddI(ixz, mim[2], cap[i], 1)
                    if nodeulv[nodenum] == '13.8':
                        ixz = getIxz(node[-2], Z[nodenum, nodenum])
                        res2[i] = decideaddI(ixz, mim[-1], cap[i], 1)

        for node in nodedata:
            nodenum = int(node[0] - 1)
            # print(node)
            for i, cnode in enumerate(capnode):
                if cnode == node[0]:
                    p = cap[i]
                    s = p / 0.98
                    q = s * 0.02
                    _q = -s * 0.02
                    r = abs(np.real(Z[nodenum, nodenum]))
                    x = abs(np.imag(Z[nodenum, nodenum]))
                    un = 1
                    du = 0.1
                    resu[i] = checkU(un, r, x, p, q, _q, du)
        for i in range(len(resu)):
            if res2[i] == True and resu[i] == True:
                res3[i] = True
        isalltrue = True
        for i in range(len(res3)):
            if res3[i] == False:
                isalltrue = False
        if isalltrue:
            isok = True
        for i in range(len(resu)):
            if res2[i] == False or resu[i] == False:
                cap[i] = cap[i] - cap[i] * 0.01
                
    #获取方式名和项目名            
    caseID,project=output_name(filename)
    
    reslist = []
    # print(res2)
    # print(resu)

    # for mkey in nodedict:
    #     hcap=cap[capnode.index(mkey)]
    #     print(hcap)

    # print(res1)
    datadict = {}

    for i, node in enumerate(capnode):
        l = {}
        l['方式名'] = caseID
        l['项目名'] = project
        l['电压等级'] = nodeulv[int(node - 1)]
        l['有功功率'] = Pi[int(node - 1), 0] * 100
        l['下级节点'] = nodedict.get(nodename[int(node - 1)])
        if cap[i] > 0:
            l['分布式光伏承载力'] = round(cap[i] * 100, 6)
        else:
            l['分布式光伏承载力'] = 0

        l['反向负载率'] = res1[i]
        l['反向负载率百分比']=str(round(res1[i], 6)*100)+'%'

        if preires[i] == True:
            l['短路电流校核'] = '通过'
        else:
            l['短路电流校核'] = '不通过'
        if preures[i] == True:
            l['电压偏差校核'] = '通过'
        else:
            l['电压偏差校核'] = '不通过'
        if preires[i] == False or preures[i] == False or round(res1[i], 2) > 0.8:
            l['评估等级'] = '红色'
        elif preires[i] == True and preures[i] == True and res1[i] <= 0.8 and round(res1[i], 2) > 0:
            l['评估等级'] = '黄色'
        elif preires[i] == True and preures[i] == True and round(res1[i], 2) <= 0.0:
            l['评估等级'] = '绿色'
        for key in nodedict:
            if nodename[int(node - 1)] in nodedict[key]:
                l['上级节点'] = key
                break
            else:
                l['上级节点'] = ''
        # print(l)
        datadict[nodename[int(node - 1)]] = l

    notred = []

    def checklv(root, mdatadict):
        """整理仅考虑节点自身情况下评估等级非红的节点

        参数
        -----------------
        root : string
            系统最高级节点
        mdatadict : dict
            节点评估结果
        """
        if mdatadict[root]['评估等级'] == '红色':
            pass
        else:
            notred.append(root)
            if mdatadict[root]['下级节点']:
                for p in mdatadict[root]['下级节点']:
                    checklv(p, mdatadict)

    checklv(rootnode, datadict)
    for key in datadict:
        if key not in notred:
            datadict[key]['评估等级'] = "红色"
    
    #问题所在
    for node in datadict:
        reason=[]
        if datadict[node]['评估等级'] == "红色":
            if datadict[node]['短路电流校核'] == '不通过':
                reason.append('短路电流校核')
            if datadict[node]['电压偏差校核'] == '不通过':
                reason.append('电压偏差校核')
            if datadict[node]['反向负载率'] > 0.8:
                reason.append('反向负载率')
            if reason == []:
                reason.append('上级节点评估等级为红色')
        if datadict[node]['评估等级'] == "黄色":
            reason.append('反向负载率')
        if datadict[node]['评估等级'] == "绿色":
            reason.append('')
        datadict[node]['问题所在'] = reason

    for key in datadict:
        # print(datadict[key]['评估等级'])
        if datadict[key]['评估等级'] == "红色":
            datadict[key]['分布式光伏承载力'] = 0.0
    for key in datadict:
        base = datadict[key]['有功功率']
        randdata = []
        for i in range(24):
            randdata.append(random.uniform(0.1, 0.2))
        randdata1 = []
        for i in range(24):
            randdata1.append((1 + randdata[i] - sum(randdata) / 24) * base)
        mysystem = platform.system()

        if mysystem == 'Windows':
            plt.rcParams['font.sans-serif'] = ['SimHei']  # win用来正常显示中文标签
        elif mysystem == 'Darwin':
            plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']  # mac用来正常显示中文标签
        plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

        plt.clf()
        x = range(24)
        y1 = np.array(randdata1)
        f1 = np.polyfit(x, y1, 3)
        p1 = np.poly1d(f1)
        yvals1 = p1(x)
        plot2 = plt.plot(x, yvals1, 'r')
        plt.xlabel('时序')
        plt.ylabel('功率/MW')
        plt.xticks(range(24))
        plt.title('节点{}日功率数据'.format(key))
        plt.savefig(global_var.get_value('direction') + 'files_generated' + os.sep + key + '.png')
        plt.clf()

    #B数据卡加上评估数据
    for i in range(0,len(bcarddata_list)):
        if bcarddata_list[i][3] in datadict.keys():
            bcarddata_list[i].append(datadict[bcarddata_list[i][3]]['分布式光伏承载力'])
            bcarddata_list[i].append(datadict[bcarddata_list[i][3]]['评估等级'])
            bcarddata_list[i].append(datadict[bcarddata_list[i][3]]['问题所在'][0])
        else:
            bcarddata_list[i].append('')
            bcarddata_list[i].append('')
            bcarddata_list[i].append('')
    # print(bcarddata_list)
    bcarddata = np.array(bcarddata_list)
    # print(datadict)
    # print(bcarddata)
    plottree(nodedict,rootnode)
    return datadict,bcarddata,lcarddata,tcarddata,nodedict,rootnode